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In the process of translation, ribosomes read the genetic code on an mRNA and assemble the 
corresponding polypeptide chain. The ribosomes perform discrete directed motion which is well 
modeled by a totally asymmetric simple exclusion process (TASEP) with open boundaries. Using 
Monte Carlo simulations and a simple mean-field theory, we discuss the effect of one or two "bot- 
tlenecks" (i.e., slow codons) on the production rate of the final protein. Confirming and extending 
previous work by Chou and Lakatos, we find that the location and spacing of the slow codons can 
affect the production rate quite dramatically. In particular, we observe a novel "edge" effect, i.e., 
an interaction of a single slow codon with the system boundary. We focus in detail on ribosome 
density profiles and provide a simple explanation for the length scale which controls the range of 
these interactions. 

PACS numbers: 05.70.Ln, 64.90.+b, 87.14.Gg 



I. INTRODUCTION 



Models and methods from nonequilibrium statistical physics find natural applications in many biological 
systems, and problems from biology have inspired many nonequilibrium models. A particularly famous 
case is translation, or protein synthesis, which motivated one of the most paradigmatic nonequilibrium 
models, namely, the totally asymmetric exclusion process (TASEP) QjSSIEIQI- In this article, we 
revisit this venerable problem and ask some new questions, related to the phenomenon of codon bias 
0, 0, IS ^1 0|' To set the scene, we give a very brief, and necessarily concise, description of how a 
protein is produced from a given gene, or to be more exactly, how the genetic information stored in 
its associated mRNA is translated into a growing polypeptide chain. Each mRNA molecule has two 
distinct ends, conventionally labelled as 3' and 5'. In the first step (initiation), a ribosome binds to the 
5' end of an mRNA. With the help of several initiation factors, the ribosome scans the mRNA until 
it encounters a start codon (usually AUG), which sets the stage for protein synthesis. This can be 
a complicated process, since the initiation efficiency depends on many factors such as growth factors, 
infections by viruses, temperature and nucleotides surrounding the start codon Then elongation 

drives translation forward, i.e., the ribosome moves codon by codon along this mRNA template until it 
reaches the stop codon, which terminates the translation process in the presence of a release factor. At 
each codon, the aminoacyl-tRNA(aa-tRNA) with the associated anticodon binds to the ribosome, and 
adds the corresponding amino acid to the growing polypeptide chain. Again, this is a complex multi-step 
process. At termination, the completed polypeptide chain is released, and the ribosome unbinds from the 
3' end of the mRNA and dissociates. Typically, at any time, several ribosomes are bound to the mRNA, 
and several protein-synthesis processes take place simultaneously. However, the ribosomes cannot overlap 
or overtake one another. The released polypeptide chain still needs to fold properly in order to function 
in a certain cell. But the focus of our research resides on the process from initiation to termination. 

To model this sequence of events in a highly simplified fashion, we start from a totally asymmetric simple 
exclusion process (TASEP), defined on a one-dimensional (ID) lattice with open boundaries, occupied by 
particles and holes. The particles jump to the right with a site-dependent rate, provided the destination 
site is empty. They enter (exit) the lattice at the left (right), with a given entrance (exit) rate. Each site on 
the lattice represents a codon on the mRNA and the particles model the ribosomes; injection, hopping, and 
drainage are associated respectively with initiation, elongation, and termination in biological terms. The 
elongation rates are commonly modeled in terms of generally accepted concentrations ("abundances") 
for the associated aa-tRNAs. While there are 64 codons, there are only about 60 anticodons 13] to 
associate with different aa-tRNAs, leading to at most 60 different hopping rates for the whole gene 
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(which typically contains hundreds to thousands of codons, i.e., sites). Moreover, there are only 20 amino 
acids, so that certain codons can be replaced by others (so-called synonymous codons) without modifying 
the final protein product. In the biological application, an important observable is the protein production 
rate, given directly by the particle current. Clearly, the protein production rate is one essential factor for 
determining the gene expression level; another one is of course simply the associated mRNA concentration. 

Clearly, this simple model falls short of the biological system in several significant aspects. One is that 
the ribosome "covers" several codons ^^j^^i as opposed to a particle occupying only a single site. A 
TASEP with extended objects was first investigated by U, and more recently in Another is that 

we model multi-step processes, such as initiation and elongation, in terms of just one rate. Therefore, we 
should not expect our findings to be fully quantitative. However, we believe studying appropriate simple 
models can provide crucial insights into universal properties that lead to useful predictions. 

In many biological or medical investigations, it is desirable to maximize or minimize the production 
of a particular protein. In the following, we focus on maximizing ("optimizing") its production, but our 
analysis can easily be applied to the opposite goal. In other words, we have to identify the rate-limiting 
step, and attempt to modify it. Here, we assume that the availability of the required aa-tRNA, rather 
than some internal reaction rate, controls the time scale on which the ribosome moves from one codon 
to the next. It is then quite intuitive, and will be shown for the model below, that the aa-tRNA with 
the lowest concentration controls the protein production rate. In order to enhance production, we can 
either over-express the rare aa-tRNA, or attempt to swap the associated codons for synonymous ones 
which employ a more abundant aa-tRNA. Here, we explore some aspects of the second mechanism; results 
for the first will be reported elsewhere pj|]. Codons associated with rare (abundant) aa-tRNAs will be 
termed "slow" ("fast"). 

Obviously, swapping all slow codons for faster ones maximizes the production rate, but will require a 
significant investment of laboratory effort. It is natural to inquire if a slightly less than maximal current 
can be attained with a much smaller amount of laboratory work. In other words, is it possible to achieve 
a significant (if not maximal) enhancement of the production rate by replacing only a small number of 
carefully targeted codons? A naive and intuitive approach would be to remove the slowest codons. But is 
the elongation rate the only factor? Or do the locations and spacings of the slow codons also play a role? 
Unless one is guided by some mathematical insight into how slow codons affect the protein production 
rate, selecting the "right codon to replace" will be a haphazard process of trial and error. 

In this article, we attempt to provide some guidance for this selection process, by considering a highly 
simplified scenario. Neglecting almost all of the inhomogeneity of the genetic sequence, we focus on a 
simple "designer gene" , consisting of many repeats of the same codon, except at one or two locations. At 
these defect sites, we insert a single slow codon. By varying the elongation rate of these special codons, 
as well as their locations and spacing, we can study their effect on the protein production rate of such a 
simple gene. 

A closely related question, namely, how to identify the rate-limiting step of the protein production 
process, was already considered in Chou and Lakatos (referred to as CL in the following) placed 
clusters of slow codons into an ordinary TASEP and varied their locations and spacings. They found that 
a single defect lowers the production rate significantly, and that a small number of slow codons, spaced 
closely together, can lower the current by an additional factor of 2 or more. The latter observation is 
interpreted as an effective interaction of slow sites with one another. Our results confirm and extend 
their findings. In particular, we present more precise data for a single slow site and discover that there is 
indeed an "edge effect", i.e., an interaction of the slow site with the system boundary, so that the particle 
current does depend on the position of the slow site. This phenomenon was not noted by CL, due to the 
larger error bars on their data. Also, we focus in more detail on ribosome density profiles and provide a 
simple explanation for the length scale which controls the range of these interactions. 

To set the stage for this investigation, we briefly review the analytical results for the steady state of 
TASEP, as well as an earlier relevant study For a TASEP with open boundaries and homogeneous 
(bulk) hopping rate 7, the steady state current is determined by the parameters a (entrance rate) and 
(3 (exit rate), both expressed in units of 7. The lattice size is denoted by N, but plays no role in the 
thermodynamic limit {N 00). The exact solutions for this model are discussed in 2, 3,0, IE 0|- Three 
phases are found: a low-density (L) phase for a < P, a < 1/2, with current a(I — a); a high-density 
(H) phase for (3 < a, (3 < 1/2, with current /3(1 — /3); and a maximal-current (M) phase for a > 1/2, 
(3 > 1/2, with current 1/4. Turning to systems with a single slow site at position k, ^3 considered a 
rather restricted case in which the defect site, with jump rate q, is located at the center of the lattice 
(fc = N/2). For large lattices {N ^ 1), an approximate stationary solution can be found by dividing the 
lattice into two separate sublattices, connected by the "defect" bond {k,k + 1). The rate q across this 
bond, along with the average occupancies at sites k and k + I, then controls the effective exit rate, /3e//, 
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from the left, and entry rate, cteff, into the right, sublattice. Focusing on the N oo limit, and using 
exact results for the usual TASEP, combined with a mean-field approximation for the current through 
the defect bond, Jg(fc), the resulting phase diagram can be determined 0. For q < 1, the phases remain 
unchanged but the phase boundaries in the phase diagram shift. One finds 

aeff = Peff = Qeff = (1) 

leading to the conditions a > /3, /? < ge// for the H phase with current — /?); I3 > a, a < q^ff for the 
L phase with current a{l — a); and finally, a, > qeff for the M phase with current 

The argument (oo) reminds us that this result is only valid if A^, fc ^ 1. We note, for completeness, that 
this simple mean-field theory can be systematically improved by considering correlations in a larger (but 
still finite) neighborhood of the slow site . 

Since we are also interested in having two defect sites (or, more precisely, bonds) in the system, it is 
natural to generalize this mean-field approach to three coupled sublattices, with the same (stationary) 
current flowing through each of them |l7j . To keep the number of parameters small, we restrict the 
discussion to a = /3 = 1. The two slow sites are placed at locations fci and ^2, separated by a distance 
d = k2 — ki, and have the same rate, q. We find that all effective exit and entry rates, for each of the 
three sublattices, are equal, again given by Qeff, Eq. As a result, the left sublattice is in an H 

phase, the right is in an L phase, and the central sublattice, characterized by aeff = (3eff, displays a 
shock, reflecting the coexistence of H and L phases. In the ordinary TASEP, such shocks are found on the 
coexistence line a — p. Their width is microscopic [T^l20l|. i.e., the density changes from the L value (a) 
to its H counterpart (1 — a) over, typically, a few lattice spacings. Moreover, they diffuse freely between 
the boundaries, so that a configurational average results in a linear density profile. Returning to our 
system with two slow bonds, these results provide us with the associated (asymptotic) current, 

provided the two slow sites are well separated. Comparing Eqs.(j2l ^^'^ ®i recognize that the second 
slow site has no further effect on the current. This statement is easily generalized to having two slow 
sites different rates qi ^ q2] in this case, the smaller rate (i.e., min{gi,g2}) sets the current through the 
system. 

To put our work in context, let us also note that a TASEP with quenched random rates on the 
entire lattice, followed by an average over the disorder distribution, was investigated by l2l|. Localized 
inhomogeneities, at the multicritical point a — (3 = 1/2, were considered by Finally, [23 extended 
the work of to extended objects. 

This paper is organized as follows. In the next section, we present our mathematical model and some 
technical details of the simulations. We then present our data and discuss the effect of defect location 
and spacing on the particle current. Finally, we turn to the implications for protein production rates and 
conclude with some simple qualitative predictions. 



II. MODEL AND METHODS 



We use a one-dimensional (ID) lattice of N sites as a template. The microconfiguration of the lattice 
can be described in terms of site occupancies^ n^, where the index i — 1, ...,N is a site label. Each site, 
initially chosen to be empty, is allowed to be occupied by a single particle (rii = 1) or left empty (n^ — 0). 
Particles enter at the left end, jump to the neighboring site on the right provided it is empty, and finally 
exit from the right end. In our random sequential updating scheme, we select a site at random and update 
it, if possible, according to the following rules: 

• ^ 1 at site 1 with rate a, 

• 1 — > at site N with rate /3, 

• 10^01 at sites (i, i -I- 1) with rate ji. 
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(a) 



(b) 



FIG. 1: Sketch of an ordinary TASEP with one slow site at fc with rate g < 1. 



For the usual TASEP, the bulk rate 7^ is chosen to be unity, for i = l,...,A'' — 1. To study how slow 
codons influence the flnal protein production rate, we modify 7^ locally, by introducing one or two slow 
sites. To introduce one slow site at position fc, we set 7^ = (7 < 1, while 7^ = 1 for i ^ k. We are 
particularly interested in the relative change in the current, Jq{k), as the location, fc, of the slow codon 
is varied. 

To introduce two slow sites, at locations fci and with separation d = {k2 — ki), we reduce both local 
rates, jki and 7^^, to g < 1. Keeping in mind that slow codons might be clustered closely together, we 
study the associated current, denoted by Jq{d), for a range of q and d. 

In our simulations, we keep a list of all occupied sites plus a single "virtual site" i = 0, which is always 
occupied and accounts for attempted particle entries into the system. To achieve the most efficient 
updating and to reduce the number of parameters in the system, we set a — l3 — ji = 1 except at one or 
two "slow" sites. At the beginning of each Monte Carlo step (MCS), we first count the number of particles 
(Af) on the lattice. Then, we choose one of the M + 1 list entries at random and attempt to update it. 
One MCS is completed after M + 1 update attempts have been made. As a result, all particles on the 
lattice and a new particle have, on average, experienced one update attempt. Typically, 5 x 10^ MCS 
are discarded to ensure that the system has reached the steady state. Results are obtained by averaging 
over least 5 x 10'* measurements, separated by 100 MCS in order to avoid correlations. Such steady state 
averages will be denoted by (...). The system size N ranges from 200 to 1000, with most data taken for 
N = 1000. 

We monitor several observables to characterize the steady state of the system. First, we measure the 
average particle current J, defined as the average number of particles entering the system per unit time. 
By the very definition of "steady state" , this current is uniform throughout the system, and could equally 
well be measured across any bond, or at the exit point. In the biological system, this current corresponds 
to the protein production rate. We also accumulated local density profiles, pi = (rii), to understand how 

they are affected by the presence of slow sites. The overall density, p = ^i'^^ follows naturally 

from these profiles. 



III. MONTE CARLO SIMULATION RESULTS 



We begin by placing one slow site (or defect bond) on the lattice as in Fig. ^ Fig. El shows several 
density profiles, illustrating the presence of significant non-uniformities. 

The "tails", i.e., the deviations from the relatively slowly varying bulk values, are quite noticeable in 
the vicinity of both the slow site and the edges of the system. Though reminiscent of the profiles shown 
schematically in ||2^, ours differ qualitatively, as a result of the loss of the i <!=;> N — i — 1 symmetry 
(fc ^ N/2), as well as a — P = 1 instead of 1/2. Not surprisingly, there is no discernable relationship 
between the profiles of the two sublattices (except in the inset). More significantly, for our case the profiles 
(within each sublattice) are non-monotonic, a feature that necessarily contradicts mean-field predictions. 
Turning to the current, we see that, except for the smallest g's, serious deviations from Eq. JSJ emerge. 
Fig. 13 for g = 0.6, shows that the current increases monotonically when the slow site is located closer 
and closer to the boundaries. Other choices of q lead to similar behavior. Since particle-hole symmetry 
holds, Jq{k) is symmetric under inversion, k ^ N + 1 — k. To quantify the /c-dependence, we define a 
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FIG. 2: Density profiles for a = 1000 lattice with one 
slow site at fe = 2 (black), 10 (dark grey, red online) and 
82 (light grey, green online) with q = 0.6. Inset: Density 
profiles for q — 0.2, 0.4, 0.6 and 0.8 (from top to bottom 
on the left, and bottom to top on the right). The slow site 
sits at the center (fc = 500), and A'' — 1000. In all cases, 
the profiles are discontinuous across the defect bond. 



FIG. 3: Jq{k) as a function of the position fcof the 
slow site for q = 0.6 and N = 1000. Jq{k) ap- 
proaches the limit 0.2463(5) as fc — > 500. The inset 
shows that Jq{k) is independent of A*', within sta- 
tistical fluctuations. 



relative change in the current, 



Ai(<z) 



J,(l) - J,(oo) 

Jq{oo) 



(4) 



As illustrated by Fig. 01 the magnitude of this difference depends sensitively on q, reaching a maximum 
at (7 = 0.49 where the relative current increase is about 2.5%. We refer to this phenomenon as the "edge 
effect". Since the current through the left and the right sublattices is controlled by the bulk densities 
there, our findings immediately imply that these bulk densities, denoted by ptuik, also shift with k. This 
feature is clearly displayed in Fig. |2J 

Returning to Fig. 13 we note that significant deviations from the limiting value, Jq{oo), are limited to a 
narrow window of (5 ~ 20 sites near the boundaries. Thanks to charge-parity (CP) invariance, both entry 
and exit edges display identical behaviors, therefore we may restrict ourselves to, e.g., the region near the 
entrance. We believe that the origin of this length scale can be traced to the presence of exponential tails 
in the density profiles of the ordinary TASEP. For a homogeneous TASEP in the H phase, with entrance 
and exit rates a and /3, the density decays exponentially into the bulk, as p( — ptuik ~ exp(--i'/^). For 
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FIG. 4: Ai(g) plotted vs. g, for iV^lOOO. 
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FIG. 5: Density profile ior q = oo, k = 500 and iV=1000. 
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FIG. 6: Density profiles for two slow sites, both 
with q = 0.2 (black), 0.5 (dark grey, red online), 
and 0.8 (light grey, green online). The system size 
is A'^ = 1000, and the two slow sites are located at 
ki = 500, and k2 = 670, resulting in d = 170. 



a > 1/2, the decay length becomes independent of a and is given by 



In [4/3(1-/3)] 



(5) 



In our case, we have a — 1, while q^ff — q/{l + q) plays the role of /3. Thus, for the q = 0.6 case, we find 
^ileff) — 15.5. If the slow site is placed so close to the boundary that k < ^{qeff), we should certainly 
expect to see deviations from Eq. (0) , a formula underpinned by the assumption k ^ 1. In support of our 
conjecture, we note that first, the observed S is consistent with £,{qeff), and second, that S, like C(9e//), 
is at most weakly dependent on the system size (cf. the inset of Fig.Oll. More detailed investigations are 
in progress, to settle this issue decisively ^17,]. In particular, if power laws such as the ones observed by 
[2^ were to prevail, this picture would have to be revised. 

According to the mean-field theory described in the former section, the presence of a defect with q > 1 
(a "fast site"), located at the center of the lattice, should have no noticeable effect on the current [Tg- Of 
course, it is not immediately apparent whether this statement remains true if the fast site is moved closer 
to the system boundaries. To explore whether such a an edge effect emerges, we consider the extreme 
case oi q = 00. Our simulation results confirms that the current does indeed remain unchanged. We find 
Jq{k) — 1/4+ 0{1/N), consistent with the expected behavior of the M phase. In contrast to the current, 
the density profiles display a dramatic signature of the fast bond, as illustrated by Fig. |5| 

If we consider the edge effect as an interaction of the slow site with the lattice boundaries, the natural 
next step is to explore the interactions between two slow sites "9^. In order to avoid edge effects, we place 
the two slow sites sufficiently far away from the boundaries and vary their separation. 

Fig. shows several typical density profiles. If q is rather small (e.g., 0.2), CL already noted the 
expected linear behavior in the central section, caused by the "wandering shock." For larger q, however, 
the center profile begins to develop distinct tails near the two slow sites. Turning to the current, Jq{d), we 
see from Fig. [Tlthat it is consistent with Eq. for d > 50, up to a finite-size correction of 0{1/N). In 
contrast, confirming the results of CL, we observe significant deviations from Eq. ((JJ, when d is decreased. 
We quantify the difference by defining 



A2(<z) 



J,(l) - J,(oo) 

Jg(oo) 



(6) 



where the arguments now refer to the distance between the two slow sites. In contrast to Ai, we observe 
that A2 exhibits a sizable dependence on q, especially for small values of q. Indeed, as already noted 
by CL, one can show that, in the limit oi q ^ the current decreases by a factor of 2. The data in 
Fig. IHlare clearly consistent with this conclusion. To sum up in words, two bottlenecks near each other 
have a dramatic effect on the current. Following CL, we may regard this phenomenon as an "interaction" 
between the two slow sites, inducing far more "resistance" when they are close than when they are 
well-separated. 
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Two additional comments are in order. First, we return to one of the predictions of the mean-field 
theory, namely that a second slow site, spaced far apart from its partner, should have no further effect 
on the current. Our data indicate that the current for two slow sites, spaced far apart, is systematically 
lower than the current for a single slow site, but only by a very small amount (less than 1%). Second, we 
can again attempt to identify a length scale which controls how Jq{d) approaches Jg(oo), as d increases. 
Since the central section of the system displays a shock, it is natural to ask whether the intrinsic width 
of the shock sets this length scale. According to 0,1131, this width covers only a few lattice spacings in 
the periodic TASEP with a single defect. Here, however, it appears that the shock broadens; preliminary 
data indicate a width of about 40 sites for q = 0.2, which is not inconsistent with Fig. [T] Again, 
this behavior appears to be independent of the system size N. More work is needed to fully explore this 
intriguing characteristic. 



IV. SUMMARY, CONCLUSIONS AND OUTLOOK 



To summarize, we study an inhomogeneous open-boundary TASEP with one or two "slow" sites {q < 1). 
Our key findings are as follows: For a single slow site, there is an "edge effect": moving the defect closer 
to the boundary enhances the current. The relative enhancement depends on q, but the effect is relatively 
small under all model conditions (at most 2.5%, for q ~ 0.49). A single fast site, on the other hand, has 
no effect on the current, irrespective of its location. A much more significant effect, with clear biological 
implications, emerges in the case of two slow sites. This was already noted by CL, and we confirm their 
findings: As a function of the separation between the two sites, the current Jq{d) decreases significantly, 
as the two sites approach each other. A quantitative measure of this effect is the fractional reduction: 
A2{q). Its dependence on q (Fig.|SJ) is nontrivial: In the g — > limit, the current is reduced by as much 
as a factor of 2. 

In order to gain a better understanding of these "interactions" between slow sites, and between slow sites 
and the system boundaries, we investigate particle (ribosome) density profiles. Every slow site displays 
a clear signature in the density profile, fully consistent with those in previous studies flslE^. If a defect 
is located at site fc, the profile is discontinuous between k and k -\- 1, and this discontinuity is surrounded 
by a "boundary layer" , or "tails" , where the densities deviate significantly from their bulk (asymptotic) 
values. In addition, the profiles display boundary layers near the system edges, as in the ordinary TASEP. 
When a defect is placed so close to a system edge or another defect that these boundary layers begin to 
overlap, the particle current develops a sensitivity to the defect-defect or defect-edge separation. In all 
other cases, the current is limited by the slowest codon in the system. In this sense, the slowest codon 
acts as a "gate keeper" . 

The above findings are significant in the sense that these currents are directly linked to the protein 
production rate. Therefore, our results should be directly applicable to "designer genes" , repeating the 
same codon, except at one or two locations. If the defect codons are "fast", i.e., associated with a highly 
abundant aa-tRNA, the production rate of the corresponding protein is insensitive to the presence of 
the defect codons, but the ribosome distribution on the mRNA will display a kink. In contrast, if the 
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FIG. 7: Jq{d) for q = 0.2 and iV = 1000, as a function of 
d. One slow site is located at ki = 500, and k2 is varied. 
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defect codons are "slow", i.e., associated with rare aa-tRNAs, the protein production rate is significantly 
reduced. The magnitude of the reduction depends on the locations of the slow codons. A single slow 
codon near the beginning (or end) of the gene allows for a higher production rate than a single slow 
codon further away, and two slow codons placed next to one another gen erate a much more drastic 
reduction than two slow codons spaced far apart. Preliminary studies |l7l | indicate that our findings 
remain qualitatively correct even if the particles (ribosomes) cover more than just one site. 

We can venture some even more wider-ranging predictions. Since several different aa-tRNA (anti- 
codons) can be associated with the same amino acid, it is possible to produce the same protein from 
several different codon sequences which will have different production rates. Given a particular gene, we 
can obviously maximize the production rate of its associated protein by systematically replacing all slow 
codons with synonymous, faster ones. However, in many genes this requires a large number of substitu- 
tions which tends to be impracticable. Instead, our findings lead us to believe that we can pinpoint a 
small number (two or three) of selected substitutions (focusing on the slowest codons, or groups of several 
slow codons clustered together) which lead to nearly optimized production rates, with considerably less 
effort. Preliminary data for real codon sequences lend first support to this conjecture, and work is in 
progress to test these ideas more thoroughly in silicon and in vitro QJj. 
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